source(here::here("code/load.R"))

dat = readRDS(here("data/data_2021-05.rds"))
setDT(dat)


##############################################
## Outcome averages and standard deviations ##
##############################################
dat[, lapply(.SD, mean), .SDcols = patterns("outcome")]
dat[, lapply(.SD, sd), .SDcols = patterns("outcome")]

# Support dummies
f <- \(x) sprintf("%s%%", mean(x[x != 5] > 5))
dat[, lapply(.SD, f), .SDcols = patterns("outcome")]


#######################
#  correlation table  #
#######################
tmp = copy(dat)[, .(coal_outcome, tax_outcome, dictator_outcome)]
setDF(tmp)
setnames(tmp, new = c("Coal", "Tax", "Amnesty"))
datasummary_correlation(
    tmp,
    escape = FALSE,
    title = "Correlations between individual-level support for buyouts across three policy domains.\\label{tab:correlation}",
    output = here("txt/fig/correlation_outcome_2021-05.tex"))


############################
## Descriptive statistics ##
############################

variable_labels = c(
    "Age" = "age_census",
    "Gender" = "gender",
    "Income" = "income",
    "Education" = "education_census",
    "Party ID" = "party",
    "Tax: Argument" = "tax_counter",
    "Tax: Recipient" = "tax_loser",
    "Tax: Support" = "tax_outcome",
    "Coal: Argument" = "coal_counter",
    "Coal: Recipient" = "coal_loser",
    "Coal: Support" = "coal_outcome",
    "Dictator: Argument" = "dictator_counter",
    "Dictator: Support" = "dictator_outcome")

tmp = copy(dat)[
    is.na(gender), gender := "Other"][
    , ..variable_labels]
setnames(tmp, old = variable_labels, new = names(variable_labels))

datasummary_skim(tmp,
    title = "Descriptive statistics for continuous variables collected in the May 2021 survey.",
    output = here("txt/fig/skim_numeric_2021-05.tex"))

datasummary_skim(tmp,
    type = "categorical",
    title = "Descriptive statistics for categorical variables collected in the May 2021 survey.",
    output = here("txt/fig/skim_categorical_2021-05.tex"))


#################
## Experiments ##
#################

tmp = copy(dat)[
    , coal_counter := factor(coal_counter)][
    , coal_counter := relevel(coal_counter, ref = "Control")][
    , tax_counter := factor(tax_counter)][
    , tax_counter := relevel(tax_counter, ref = "Control")][
    , dictator_counter := factor(dictator_counter)][
    , dictator_counter := relevel(dictator_counter, ref = "Control")][
    , coal_worker := fcase(
        coal_loser == "Workers", 1,
        coal_loser == "Corporations", 0)][
    , tax_worker := fcase(
        tax_loser == "Workers", 1,
        tax_loser == "Corporations", 0)]


cr = c(
    "dictator_counterGovernment intervention" = "Government intervention",
    "dictator_counterMoral hazard" = "Moral hazard",
    "dictator_counterMoral principle" = "Moral principle",
    "coal_counterGovernment intervention" = "Government intervention",
    "coal_counterMoral hazard" = "Moral hazard",
    "coal_loserWorkers" = "Workers",
    "tax_counterGovernment intervention" = "Government intervention",
    "tax_counterMoral hazard" = "Moral hazard",
    "tax_loserWorkers" = "Workers",
    "coal_counter" = "Counter",
    "coal_loser" = "Target",
    "tax_counter" = "Counter",
    "tax_loser" = "Target",
    "dictator_counter" = "Counter"
)

mod = list(
    "Coal" = lm(coal_outcome ~ coal_counter + coal_loser, data = tmp),
    "Tax" = lm(tax_outcome ~ tax_counter + tax_loser, data = tmp),
    "Dictator" = lm(dictator_outcome ~ dictator_counter, data = tmp))

# reference requested by R1
coef(mod[[1]])["coal_loserWorkers"] / mean(dat$coal_outcome[dat$coal_loser == "Corporations"], na.rm = TRUE)

modelsummary(mod,
    output = here("txt/fig/losers_counters_2021-05.tex"),
    title = "Estimated average treatment effects in the Coal, Tax, and Dictator vignettes.\\label{tab:losers_counters}",
    escape = FALSE,
    vcov = "HC1",
    align = "lccc",
    coef_rename = cr)

bkg = list(geom_vline(xintercept = 0, linetype = 3))
p = modelplot(mod,
    vcov = "HC1",
    coef_map = cr,
    color = "black",
    background = bkg) +
    facet_grid(~model) +
    xlab("Estimates and 95% intervals") +
    theme_vab()

# ggsave(p, file = here("txt/fig/losers_counters.png"), width = 5, height = 2)
ggsave(p, file = here("txt/fig/losers_counters_2021-05.pdf"), width = 5, height = 2)



# Replication with Lin (2013) interaction terms
mod = list(
    "Coal" = lm(coal_outcome ~ (coal_counter + coal_loser) * (age_census + education_census + gender), data = tmp),
    "Tax" = lm(tax_outcome ~ (tax_counter + tax_loser) * (age_census + education_census + gender), data = tmp),
    "Dictator" = lm(dictator_outcome ~ (dictator_counter) * (age_census + education_census + gender), data = tmp))

mfx <- list()
mfx[["Coal"]] <- avg_comparisons(mod[[1]], variables = c("coal_counter", "coal_loser"), vcov = "HC3")
mfx[["Tax"]] <- avg_comparisons(mod[[2]], variables = c("tax_counter", "tax_loser"), vcov = "HC3")
mfx[["Dictator"]] <- avg_comparisons(mod[[3]], variables = "dictator_counter", vcov = "HC3")
mfx[[1]]$term <- mfx[[1]]$contrast
mfx[[2]]$term <- mfx[[2]]$contrast
mfx[[3]]$term <- mfx[[3]]$contrast
mfx[[1]]$Issue <- "Coal"
mfx[[2]]$Issue <- "Tax"
mfx[[3]]$Issue <- "Dictator"
mfx <- rbindlist(mfx)

p <- ggplot(mfx, aes(x = estimate, xmin = conf.low, xmax = conf.high, y = contrast)) +
    geom_pointrange() +
    geom_vline(xintercept = 0, linetype = 3) +
    facet_wrap(~Issue) +
    labs(y = "", x = "Estimates and 95% intervals.") +
    theme_vab()
ggsave(p, file = here("txt/fig/losers_counters_lin_2021-05.pdf"), width = 5, height = 2)



################################################
## Socio-demographics and support for buyouts ##
################################################

tmp = copy(dat)[
    , Education := education_census][
    , Age := age_census][
    , Party := party][
    , Woman := fcase(
        gender == "Woman", 1,
        gender == "Man", 0)]

mod = list(
    "Coal" = lm(coal_outcome ~ Age + Education + Woman + Party, tmp),
    "Tax" = lm(tax_outcome ~ Age + Education + Woman + Party, tmp),
    "Dictator" = lm(dictator_outcome ~ Age + Education + Woman + Party, tmp))

cr <- \(x) gsub("Party|Education|Age", "", x)

bkg = list(geom_vline(xintercept = 0, linetype = 3))
p = modelplot(mod,
    vcov = "HC1",
    coef_omit = "Intercept",
    coef_rename = cr,
    background = bkg) +
    scale_color_manual(values = tol_palette) +
    guides(color = guide_legend(reverse = TRUE)) +
    xlab("Estimates and 95% intervals") +
    theme_vab() +
    theme(legend.title = element_blank())

# ggsave(p, file = here("txt/fig/socio_demo_determinants_2021-05.png"), width = 5, height = 2)
ggsave(p, file = here("txt/fig/socio_demo_determinants_2021-05.pdf"), width = 5, height = 3)



